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The mechanism of phase synchronization between uncoupled limit-cycle oscillators induced by 
common external impulsive forcing is analyzed. By reducing the dynamics of the oscillator to a 
random phase map, it is shown that phase synchronization generally occurs when the oscillator is 
driven by weak external impulses in the limit of large inter-impulse intervals. The case where the 
inter-impulse intervals are finite is also analyzed perturbatively for small impulse intensity. For 
weak Poissonian impulses, it is shown that the phase synchronization persists up to the first order 
approximation. 

I. INTRODUCTION 

When a limit-cycle oscillator is driven weakly by the same temporal sequence of a fluctuating input, its phase tends 
to exhibit the same dynamics repetitively among different realizations even if small external disturbances distinguish 
each realization. For example, a cortical neuron generates spikes more reproducibly when it receives a fluctuating input 
current rather than a constant input current QLEU^. This fluctuation-induced reproducibility of a single oscillator can 
be interpreted as phase synchronization between uncoupled oscillators driven by common external forcing, because 
repeated measurements on a single oscillator using the same input is equivalent to a single measurement on multiple 
identical oscillators. It indicates the existence of a physical mechanism that statistically stabilizes the limit-cycle orbit 
in the phase direction by a fluctuating input. 

There have been a variety of studies regarding this phenomenon 0, H, IE 01 H 01 EH EH Ell ■ Among them, Teramae 
and Tanaka |lfjj | made significant progress in understanding its universality from the viewpoint of nonlinear dynamics. 
Using the Stratonovich-Langevin equation resulting from the phase reduction method 0, IfH Il5j| , they generally 
proved that limit-cycle oscillators always synchronize in phase when they are driven by a vanishingly weak Gaussian- 
white forcing (see also Goldobin and Pikovsky |12|). Independently, we also analyzed this phenomenon in a different 
setting ^]|. We assumed a simple random telegraphic forcing to the oscillator that switches between two values 
randomly, which is not necessarily vanishingly small. We reduced the dynamics of the system to a pair of random 
maps, and generally showed that the oscillators always synchronize in phase when the phase map is monotonic. 

In this paper, we consider yet another model of fluctuation-induced phase synchronization. Specifically, we assume 
the external forcing to be random impulses. Such a model has wide applications to various natural phenomena, since 
impulsive noises are abundant in nature |l6j . For example, a cortical neuron interacts with other neurons via spike 
trains, which are modeled as impulses in the simplest approximation |24| . Within this model, we can generally prove 
that the oscillators actually undergo fluctuation-induced phase synchronization in the limit of large inter-impulse 
intervals. In addition, we can also discuss the case where the inter-impulse interval is finite within this model. Both of 
the previous analyses in Refs. [Tol 1 1 1 1 ] assumed that the phase distribution of the oscillator is completely uniform on 
the limit cycle, which corresponds to the assumption of vanishingly weak or infinitely slow-switching external forcing. 
However, in practical situations, such assumptions may not be valid, and the phase distribution would generally be 
non-uniform. Thus, the effect of non-uniform phase distribution on the phase synchronization should be assessed. By 
developing a perturbation theory for weak impulse intensity, we discuss the effect of slight non-uniformity of the phase 
distribution on the phase synchronization. Especially, we will show that the phase synchronization persists even if 
the phase distribution becomes slightly non-uniform for the Poissonian impulses. 

Our analysis presented in this paper will extend the class of fluctuation-driven limit-cycle oscillators that exhibit 
phase synchronization, and will provide deeper insight into this phenomenon. 
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This paper is organized as follows. In section II, a general model of impulse-driven limit-cycle oscillators is intro- 
duced, and phase synchronization is demonstrated using two typical limit-cycle oscillators. In section III, reduction of 
the dynamics of impulse-driven oscillators to a random phase map is presented. In section IV, stability in the phase 
direction is analyzed in the case where the phase distribution is uniform. In section V, effect of non-uniformity of the 
phase distribution is analyzed perturbatively for small impulse intensity. Finally, Section VI gives the summary. 
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FIG. 1: Synchronization of 50 Stuart-Landau oscillators driven by external impulses, (a) Zero-crossing events at a = and (b) 
at a = 0.1. (c) Time sequence of the averaged modulus (R) of the order parameter calculated at a — 0, 0.05, 0.1 and 0.2. 



II. PHASE SYNCHRONIZATION INDUCED BY EXTERNAL IMPULSES 



We first present a general model of limit-cycle oscillators driven by a common external impulsive forcing, which we 
will analyze in later sections. Then, before going into a general theory, we numerically demonstrate phase synchro- 
nization induced by external impulses using two typical models of limit-cycle oscillators, and briefly comment on its 
mechanism. 



A. General Model 



We consider an ensemble of N identical limit-cycle oscillators subject to a common external impulsive forcing in 
the following general form: 



X i (t)=F(X i )+I(t) 



(1) 
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for i = 1, • • • , N, where Xj(t) represents the internal state of the i-th oscillator at time t and F its dynamics. We 
assume that Eq.Q) has a single stable limit-cycle solution X (i) in the absence of external impulsive forcing I(t), 
whose basin of attraction is the entire phase space except some unstable fixed points. External impulsive forcing I(t) 
is given by 

oo 

i(t) = ^^(t-g, (2) 

n=l 

where t\, i 2 , • • ■ are occurrence times of the impulses, and e 1; e 2 , • • • are random vectors representing the "direction" 
of the impulses. When the oscillator receives an impulse at t — t n , its state Xj is suddenly kicked by a random 
displacement e n to a new state X^ + e„. 

We assume that the direction of the impulse e is mutually independent and identically distributed. We denote 
its probability density by Q(e), which is normalized as J deQ(e) = 1. We also assume the interval T between two 
successive impulses to be independent and identically distributed. We denote its probability density function by P(T), 
which is normalized as J °° dTP(T) — 1. We further assume the intervals to be sufficiently long, so that the orbit 
kicked away from the limit cycle by an impulse can return to the limit cycle before the next impulse. This is the 
necessary condition for the phase description of the oscillator, which we adopt in the following discussion. The time 
needed for this process of course depends on the characteristic of the oscillator and on the intensity of the impulses, 
which is very roughly of the order of the period of the limit cycle. 

In this paper, we mainly consider impulses generated by a Poissonian process. Then Eqs.Q) and (J2J describe a 
Poisson-driven Markov process \l& In the Poissonian process, an impulse is generated with probability v in an 
infinitesimal time interval At. The probability density P(T) of the inter-impulse interval T is given by the exponential 
distribution 

P(T) = ±exp(-^), (3) 

where t = v^ 1 determines the mean inter- impulse interval. Of course, in this Poissonian case, there exists a certain 
probability that the generated interval becomes shorter than the above-mentioned return time of a kicked orbit to the 
limit cycle. In such a case, the phase description fails to approximate the true dynamics precisely. However, when 
r is sufficiently large, such probability becomes small, and the phase description would be a good approximation 
statistically. 

The probability density of the impulse direction Q(e) should be chosen properly depending on the problem un- 
der consideration. For example, when we consider neural oscillators, usually only the membrane potential can be 
stimulated experimentally by a current injection. Therefore, the random vector e has only one non-zero element cor- 
responding to the voltage component of X, and we only need to consider one-dimensional probability density Q(cr), 
where a is the intensity of the current impulse. In the following examples, we only treat the cases where the stimulus 
can take a single fixed direction and intensity eo, namely Q(e) = S(e — eo), but we present our theory in a general 
form so that it is also applicable to the case where the stimulus takes various direction and intensity. 



B. Examples 

1. Stuart-Landau oscillator 

Our first example is an ensemble of noisy Stuart-Landau oscillators |l4j driven by a common sequence of Poissonian 
random impulses with fixed intensity. The Stuart-Landau oscillator is the simplest limit-cycle oscillator, which is 
derived as a normal form of the super-critical Hopf bifurcation jl4j . The model is described by 

Ci(t) = (l + «*)a--(l+*c 2 )|Q| 2 Q + + (4) 

for i = 1, • • • , N, where Ci is the complex amplitude of the i-th oscillator, Co and c 2 are real parameters, I(t) represents 
a sequence of external impulses, and Q(t) is a mutually independent complex Gaussian- white noise additionally 
introduced to represent small external disturbances. For simplicity, we drive only the real part of Ci by the external 
impulsive forcing I(t), which is given by 



oo 

I(t)=aY,t(t-t n ), (5) 

n=l 
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where the real parameter a represents the intensity of the impulses, t n the time of occurrence of the impulse, and 
8 the Dirac's delta function. The impulses are generated by a Poisson process with mean inter-impulse interval r. 
When the oscillator receives an impulse, its real component ReC; suddenly jumps by an amount a. It is assumed that 
the complex Gaussian- white noise Q(t) has zero- mean, and whose variance is given by 

(ReCi(t)ReG(O) = D5 UJ S(t - t'), 
(ImCi(*)ImO(0> = DS itj S(t-t'), 

(ReC»(*)ImO(0> = 0, (6) 

where D determines the noise intensity. We fix the parameters at Cq — 2, c 2 = — 1, r = 2, and the noise strength at 
Vd = 10~ 3 . 

We initially set the phase 9i of each oscillator uniformly and randomly on [0, 1] (see the next section for the precise 
definition of the "phase" ) , where the zero-crossing point of Cj from ImQ < to ImQ > is chosen as the origin of 
phase, &; L — 0. We then evolve the oscillators under the influence of the Poissonian impulses and the weak Gaussian- 
white noise. Figurcsnja) and (b) plot zero-crossing events of Ci of N = 50 Stuart-Landau oscillators by bars (so-called 
raster plot) after transient for the cases a — and a = 0.1. It can be seen that the oscillators synchronize in phase 
when a = 0.1, whereas they do not synchronize at all when a = 0. To quantify the degree of synchronization, we 
introduce an order parameter |l4j 

1 N 

i?exp(2™6) = -jyX; 

exp(27ri6>i) (7) 

»=i 

using phase &i of each oscillator. The modulus R of this order parameter takes R = 1 for complete synchronization 
and R = for complete desynchronization. Figure^c) displays temporal evolution of the modulus (R) averaged over 
50 realizations from different initial conditions. It gradually increases from a small value to 1 when a = 0.05 and 
a = 0.1, while it constantly takes a small value when a = 0. Thus, the uncoupled Stuart-Landau oscillators driven 
by a common sequence of Poissonian impulses synchronize in phase even if small external disturbances exist. 



2. Hodgkin- Huxley model 

Our second example is an ensemble of the Hodgkin-Huxley neural oscillators |24j driven by a common sequence of 
Poissonian random impulses with fixed intensity. It is given by the following set of equations |24| : 

C m Vi(t) = G Na mlhi(E Na - Vi) + G K ni(E K - V t ) + G m (V rest - V t ) + I + I(t) + 

rhi(t) = a m (l - rut) - f3 m n%i, 

hi(t) = a h (l - hi) ~ p h hi, 

rii[t) = a n (l - m) - n m, (8) 

for i = 1, • • • , TV, where Vi represents the membrane potential of the z-th neural oscillator, rrij and hi the activation of 
its sodium channel, and n.; the activation of the potassium channel, Iq the constant input current, I(t) the external 
impulsive forcing, and the additional external disturbances. Parameters Gjva, Gk, and G m represent conduc- 
tances of the channels, E^ ai and Ex represent their reversal potentials, and V rest represents the rest voltage. a x and 
/3 X (x = m, h, n) are rate constants that are given by the following equations: 

= ex P( _ Yg)' 
i 



0.1(25- 


v) 


Pm 


exp(^) 


-1' 


0.07exp(- 


--) 
20 j ' 


Ph 


0.01(10- 


■v) 


Pn 


exp(i^) 


-1' 



0.125exp(~^). (9) 

The parameters are fixed at the standard values presented in the textbook |2^|, i.e., GNa = 120, Ej^ a — 115, Gk — 36, 
Ek = —12, G m = 0.3, V res t = 10.613 and C m = 1. We fix the constant input at Iq = 11. When the external impulsive 
forcing I(t) is absent, this model exhibits stable limit-cycle oscillation. The external impulsive forcing is given by 

oo 

I(t)=aJ2 S ( t - t n), (10) 

n=l 
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where a determines its intensity. The impulses are generated by a Poissonian process with mean inter-impulse interval 
r = 100. The external disturbance is represented by a Gaussian- white noise of mean and variance D. We fix 
the noise strength at y/~D = 10~ 3 hereafter. We define the zero-crossing event ("firing event") of this Hodgkin-Huxley 
neural oscillator as the moment at which the variable V^changes its sign from Vi < to Vi > 0. We take this point 
as the origin, and define a phase along the limit cycle [ljj, [lj] . 

As in the previous case, we set the initial phases of the Hodgkin-Huxley oscillators randomly, and evolve them 
under the effect of Poissonian impulses and Gaussian- white noises. Figures |2fa),(b) show zero-crossing events of an 
ensemble of 50 Hodgkin-Huxley neural oscillators without or with external impulses (a — or a — 2) after the initial 
transient. The zero-crossing event well coincides with each other when a = 2, i.e. the oscillators synchronize in phase 
after the initial transient. Of course, they do not synchronize at all when a = 0. Figurc[2Ic) shows temporal evolution 
of the modulus (R) of the order parameter averaged over 20 independent realizations at several values of the impulse 
intensity a. When a takes 1, 2, or 4, phase synchronization occurs and (R) gradually increases from a small value to 
1. 

It is also possible to observe impulse-induced desynchronization by choosing the impulse intensity appropriately as 
shown in Fig.^a), where the zero-crossing events of 50 Hodgkin-Huxley neural oscillators are plotted. In this case, the 
external constant input is set at Iq = 9.8, the intensity of the external impulse at a = 6, and the inter-impulse interval 
at r = 50. The phase of each oscillator is initially set at roughly the same value, except tiny additional fluctuations of 
order 0.1. To avoid spurious complete synchronization due to numerical cutoff, small external Gaussian- white noise 
with \f~D = 10~ 3 is additionally applied. The phases of the oscillators scatter occasionally, and correspondingly the 
order parameter drops as shown in Fig. Efb). Thus, initial tiny phase differences among the oscillators can also be 
enhanced by the external impulses. 



(a) 

50 
| 40 
o 30 

a 20 

O 10 






y [i ,f ,/ ip* ,/ 


•f <f •! •{ •? 'i' 


,J' .]/ ' ,V ,;' _ 


-• , ' '• , ' '• 


1 ,' 1 ,' 1 .' 1 ,' 1 ,' 1 ,' 
' '• ,' '• , ' '• ,' '• , ' '• ,' '• 
'• , '• , '• , '• . '•' , '•' , 


1 ,' 1 ,' 1 ,' 1 ,' 1 ,' 1 ,' 1 ,' 
' '• ,' '• ,' '• .' '• , ' '• ,' '• ,' '• 
• , '• , '• ' , '• , '• . , , 


1 ,' 1 ,' 1 ,' 1 ,' 1 ,' 

• 


T (\ >!' 


'•'] '-'t *\ '•''] '-\ ''\ 


'■"'i ^'i ( 'i *\ > '-\ '''i 


K C\ (\ '.'\ ? 











39700 



39800 



39900 



40000 



(b) 

50 



Z^O 
2 30 

= 20 

u 

% 10 



I I I I 1 I I I II I I I 1 I I I II I I I 1 I I I 

j || j | ;i S ,i I ' 



o 

39700 



(c) 



Jl 



39800 



39900 



40000 




10000 



20000 30000 40000 

t 



FIG. 2: Synchronization of 50 Hodgkin-Huxley neural oscillators driven by external impulses, (a) Zero-crossing events at a = 
and (b) at a = 2. (c) Time sequence of the averaged modulus (R) of the order parameter calculated at a = 0, 1, 2 and 4. 
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FIG. 3: Desynchronization of 50 Hodgkin-Huxley neural oscillators induced by external impulses, (a) Zero-crossing events, (b) 
Temporal sequence of the modulus R of the corresponding order parameter. 



C. Some comments on the mechanism of impulse-induced phase synchronization 

The mechanism of phase synchronization induced by common external input is basically a single-oscillator problem, 
though we consider an ensemble of oscillators in the above examples. The origin of the phase synchronization is 
the local stabilization of each limit cycle in the phase direction due to the external impulses. Namely, small phase 
disturbances of a single oscillator shrink statistically, as we formulate in the following sections by reducing the 
dynamics of each oscillator to a random phase map. At the same time, it indicates the suppression of small difference 
in phase between any pair of oscillators. Due to random impulses, the phase of each oscillator diffuses on the limit 
cycle in addition to the constant rotation. Once two phases of any pair of the oscillators come close accidentally, 
their difference can no longer grow but shrinks statistically due to the local stability, leading to the synchrony of the 
entire ensemble. This mechanism has certain similarity to that of chaos synchronization induced by common random 

forcing [ii is mm. 

In the second example, we demonstrated that external impulses do not only lead to phase synchronization but can 
also cause phase desynchronization. Though we mainly focus on phase synchronization in this paper, this fact is 
important in understanding that fluctuation-induced phase synchronization is not a trivial phenomenon but has some 
subtleties. 



III. REDUCTION TO A RANDOM PHASE MAP 



In order to analyze the stability against phase disturbances, we first reduce the dynamics of our impulse-driven 
limit-cycle oscillator to a random phase map. 



A. Random phase map 

Following the standard procedure 0,0], we define a phase 8 — 6(Xq) along the limit cycle orbit Xo(£) in such a 
way that 9 increases with a constant angular velocity oj. The phase is normalized by the period of the limit cycle, 
so that its range is [0, 1] where and 1 represent the same phase. This definition of phase can be extended to the 
entire phase space except at phase-singular points, yielding a phase field #(X). It is achieved by identifying a point 
P in the phase space with a point Q right on the limit cycle in such a way that the two orbits started from P and 
Q asymptotically coincide. A set of points that have equal phase is called an isochron. The entire phase space is 
composed of isochrons with various phases. 
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In the absence of external impulses, the phase 9(t) = 8(X.(t)) obeys 

e(t) = u (11) 

on the entire phase space (except at phase-singular points). When the external impulses are given, the orbit is 
perturbed. Let us assume that the orbit is on the limit-cycle at time t = t n — , i.e. immediately before the n-th 
impulse (we say the orbit is "on" the limit cycle when it is sufficiently close to the limit cycle.) We denote its location 
by Xo(i„ — 0) and its phase by 9 n = 9(t n — 0). We also denote the interval between the n-th impulse and the next 
(n + l)-th impulse by T n = t n+ i — t n . When the oscillator receives an impulse e n S(t — t n ) at t = t n , the orbit is kicked 
off the limit cycle and jumps to a new phase-space point as 

X(t„ + 0) = X Q (t n - 0) + e n . (12) 

This new phase-space point X(t„ + 0) is on a certain isochron of the limit cycle, whose phase we denote by <f) n (unless 
it is kicked exactly onto the phase-singular point, which rarely occurs). We represent this mapping from 9 n to </>„ by 

n =F(0 n ,e„) = 9 n + G{9 n ,e n ), (13) 

which we call a "phase map" hereafter. In the second expression, we split F(8 n ,e n ) into the trivial part 9 n that 
exists even without any impulses, and the non-trivial part G(9 n , e„) arising from the impulse. Since F(9, e) is a phase 
map, it is a periodic function on [0,1]. Therefore, F(l,e) = F(0, e) and G(0,e) = G(l,e) should hold (we should 
treat them in modulo 1). As we discuss later, the above rule gives rise to an impulse-driven phase equation of the 
Ito-type 0[l3. 

After the arrival of the n-th impulse, the oscillator evolves freely with no external impulses from t = t n + to 
t = t n +i — for an interval of T n = t n +i — t n , and the phase changes from </>„ to <f) n + ioT n during this interval. If T n 
is sufficiently large, the orbit evolves from X(i n + 0) to a new point Xo(t n +i — 0) on the limit cycle. 

Thus, corresponding to the evolution of the variable Xo(t„ — 0) — > X(t„ + 0) — > Xo(t n +i — 0) from t = t n — to 
t = fn+l — 0, the phase evolves as 9 n — > <f> n — > <p n + u>T n . Hence we obtain the following evolution equation of the 
phase: 

e n+1 = uT n + F(6 n ,e n ) = 9 n + cuT n + G(6 n ,e n ). (14) 

Since T n and e n are random variables whose probability density functions are given by P(T) and Q{e) respectively, 
this equation describes a random map. When we consider Poissonian random impulses, the time step n roughly 
corresponds to the real time t as n ~ t/r, because the mean inter-impulse interval is r. 
If we go back to the continuous description, the dynamics of the phase 9 can be written as 

oo 

d(t)=u + J2G(6 n ,e n )5(t-t n ). (15) 

n=l 

The external impulse is now explicitly multiplicative in this equation. This impulse-driven phase equation is of Ito 
type |lfil Il7l |. namely, G{9 n ,e n ) depends only on the phase 9 n before the n-th impulse, which stems from the rule we 
have assumed for the phase jump caused by an impulse. 



B. Relation to the phase response function 

According to the standard theory of phase reduction [l3L ll4T| , when the orbit on the limit-cycle X at phase 9 
is kicked by a weak impulsive force p to another isochron, its new phase 4> is given by a linear projection of the 
perturbation p on the gradient of the phase field 0(X) as 

c/) = 9 + Z(9) ■ p, (16) 

where 

Z(0) = Vx0(X)| XoW (17) 

is the conventional phase response function representing the gradient of 0(X) on the limit cycle orbit Xo(f). Comparing 
this equation with Eqs.lO and ifHjl. we obtain 

F(6 n ,e n ) = n + Z(0 n ) - e„, G(6 n , e n ) = Z(0 B ) • e n , (18) 

so that the phase dynamics can be described by 

9 n +i = 9 n + uT n + Z(0„) • e„. (19) 

Thus, for sufficiently small e n , the phase map can simply be represented using the inner product of the conventional 
phase response function Z(0) and the direction of the impulse e„. 
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C. Generalized Frobenius-Perron equation 

Temporal evolution of the probability density function (PDF) p(0, n) of the phase at time step n is described 
by a generalized Frobenius-Perron equation 18], which is convoluted with a transition kernel W(9) that represents 
random shifting on the limit cycle for a random duration T drawn from P(T), and is also averaged by the probability 
density Q(e) of impulse directions e, 

p{8,n + 1) = f d4>W{B - 4>) f deQ{e) f d#fa - F(ip, e))p(ip, n), (20) 
Jo J Jo 

where the argument 9 — <f> of W{0 — <p) should be interpreted in modulo 1. In deriving this equation, we utilized the 
fact that 6 n , e„, and T n are mutually independent (0 n depends only on ei, ■ • • , e„_i and T\, ■ ■ ■ , T„_i) [l8l | . 

For Poissonian impulses, the explicit form of the transition kernel W(0) can be calculated from Eq.@ by taking 
into account the periodicity in and the Jacobian of the transformation, which is given by 

1 °° ft 4- -i pS/i^T) A 

w(o) = -Y, p (—^—n tj^\ = t JL -a^ m (o^ 1 )' ( 21 ) 

3=0 y ' 

where we defined A = 1/wr. Of course, it is normalized as L d0W{8) = 1. 

Sufficiently after the initial transient, the PDF p(0, n) is expected to reach a stationary state p(9), but it is generally 
difficult to calculate this stationary PDF analytically even if the map F(0, e) has a simple functional form. In the 
following, we first analyze the limit of large inter-impulse interval r where the stationary PDF becomes uniform, and 
then analyze the deviation of the stationary PDF from the uniform density perturbatively for small impulse intensity. 

D. Examples of phase maps 

1. Stuart-Landau oscillator 

Figure Ufa) plots numerically calculated phase maps F(0,<j) of the Stuart-Landau oscillator at several values of 
the impulse intensity a. As \a\ becomes larger, the phase map deforms from the trivial identity map noticeably, and 
finally becomes non-monotonic when a ~ ±0.8. Figure Iffib) displays the phase response normalized by the impulse 
intensity [F(9,a) — 0\ja = G{9,o)/o at several values of a. If a is sufficiently small, Eci. (|18fl should hold, and the 
different curves corresponding to different values of a should collapse. The limiting curve at a — > gives the real 
component of the phase response function Z(6*). It can be analytically calculated for the Stuart-Landau oscillator as 
Zj^ e (0) = sin(27r6' + 37r/4)/3-\/2 [li|. which is also shown in the figure. 

2. Hodgkin-Huxley model 

FigureEfa) shows numerically calculated phase maps F(0, a) of a Hodgkin-Huxley neural oscillator at several values 
of the impulse intensity a. As \a\ becomes larger, the phase map deforms from the trivial identity map and finally 
becomes non-monotonic at a ~ ±5. Figure |SJb) shows normalized phase response [F(9, a) — 0] /a = G(0,a)/a at 
several values of a. It can be seen that the curves actually collapse at small a, and deviates at larger a. The liming 
curve at cr — ► gives the ^-component of the phase response function Zy(6). The curve corresponding to a = 0.1 in 
Fig. GJb) gives an approximation to the phase response function. 

IV. STABILITY IN THE PHASE DIRECTION 

Synchrony of uncoupled oscillators induced by external impulses is the result of statistical stabilization of each 
oscillator against phase disturbances. Such stability is characterized by the Lyapunov exponent of the random phase 
map Eq. (jH|> . 
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FIG. 4: (a)Phase maps F(9, a) of the Stuart-Landau oscillator at several values of the impulse intensity a. (b) Normalized 
phase response [F(9) — 9] /a at several values of the impulse intensity a. Theoretical phase response function Z-^ e (9) = 
siri(27r# + 37r/4)/3\/2 that corresponds to the limit of a — > is also shown. 




FIG. 5: (a) Phase map F(8, a) of the Hodgkin-Huxley neural oscillator at several values of the impulse intensity a. (b) 
Normalized phase response [F(0) — 9] ja at several values of the impulse intensity a. The curve corresponding to a = 0.1 gives 
an approximation to the phase response function Zy{&) of the Hodgkin-Huxley neural oscillator. 



A. Lyapunov exponent 



Let us consider the temporal evolution of a small deviation A9 n from the original orbit 9 n . The linearized evolution 
equation of this small deviation is given by 



M n+1 = F'(d n ,e n )A9 n , 
where F'(0 n ,e) = (dF(9,e)/dd) g=e . At large time step n, A9 n expands as 

A9„ 



(22) 



A9 n 



H \F'{9 m ,e m )\ 



m=0 
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. xj) ^2 log\F'(e m ,e m )\ 

_m=0 

~ exp(An), (23) 
where we introduced the Lyapunov exponent A of the map F(0, e) averaged over the PDFs Q(e) and p(6), 



\ = {{log\F'(8,e)\)) e ,, 



ddp(6) / deQ (e) log \F' (9, e) 



(24) 



If A is negative, A9 n shrinks on average, so that the deviation from the original orbit is suppressed, whereas if A is 
positive, small external disturbances will be enhanced. Thus, the value of A gives a (local) condition for the phase 
synchronization. 



B. Limit of large inter-impulse intervals 

As we mentioned previously, it is difficult to obtain the stationary PDF p(9) analytically. However, when the 
inter-impulse interval r is sufficiently large, it can be approximated by a uniform density. In the limit of large r, 
the transition probability tends to be uniform, i.e. W{9) — > 1, which can easily be confirmed from Eq. (|21|l in the 
Poissonian case. Correspondingly, the stationary phase PDF p{9) approaches a uniform density in the large r limit: 



P(0) 



1. 



(25) 



In this limit, we can obtain a sufficient condition of phase synchronization for general limit-cycle oscillators: when 
the phase map F(9, e) is a monotonically increasing function of 9, the Lyapunov exponent A is always non-positive, 
namely, when F'(9, e) = 1 + G'(9, e) > is satisfied. We can then bound A from above as 



A = 



< 



d9p(9) / deQ (e) log F' (9, e) 



deQ(e)F'{6,e) 



d9p(9) log 



d9p{9) log 



J deQ{e)G'(9,e) 



< 



f d0p(6) J deQ{e)G'{8,e). 



(26) 



In the above transformation, we utilized Jensen's inequality J deQ(e)log[g(e)] < log[ J deQ (e)p(e)] that holds for a 
concave function log[- • •], the normalized probability density Q(e), and a positive scalar function g(e). The second 
inequality follows from log(l + x) < x. By using the facts that p(9) = 1 and G(0,e) = G(l,e), the upper bound of A 
can be calculated as 



d9p{9) / deQ{e)G'(6,e) 



deQ(e) f d9G'(9,e) 
Jo 

J deQ(e) [G(l, e) - G(0, e)] = 0. 



(27) 



Thus, for monotonically increasing F(9,e), the Lyapunov exponent A is always non-positive. The equality A = 
holds only when F(9,e) is a trivial identity map for all e, i.e. F(9,e) = 9, which follows from the equality condition 
of Jensen's inequality. Therefore, small deviations from the original orbit always shrink by applying random external 
impulses with large inter-impulse intervals, when the phase map F{9, e) is monotonic. 

As we mentioned previously, when e is small, G(9,e) can be represented using the phase response function Z(9) as 
G(9,e) ~ Z(#) • e. Since F(9,e) = 9 + G(9,e), F(9,e) is monotonically increasing with respect to 9 for sufficiently 
small e. Therefore, when the intensity of external impulses is small and the mean interval between impulses is large, 
A always becomes negative. 
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C. Examples 

1. Stuart-Landau oscillator 

As can be seen from Fig.0Ja), the phase map F(9,a) of the Stuart-Landau oscillator is monotonic as long as \a\ is 
small. Therefore, the Stuart-Landau oscillator exhibits phase synchronization induced by external impulses at such 
values of \a\ for sufficiently large inter-impulse intervals, as we demonstrated previously. 

2. Hodgkin-Huxley model 

Similarly, as shown in Fig.J^a), numerically calculated phase maps F(9,a) of the Hodgkin-Huxley neural oscillator 
are monotonic when |er| is not so large. Therefore, the Hodgkin-Huxley neural oscillators also exhibit impulse-induced 
phase synchronization for such values of \cr\. When |<r| becomes large, the phase map can become quite complex, 
which can lead to the impulse-induced phase desynchronization mentioned previously. 

V. EFFECT OF NON-UNIFORM PHASE DISTRIBUTION 

In the previous section, we discussed the limiting case of large inter-impulse intervals, where the stationary PDF 
of the phase becomes uniform. If the mean inter-impulse interval is not so large, the stationary PDF would generally 
become non-uniform. In this section, we first develop a perturbation theory to approximate the non-uniform PDF 
for weak external impulses. We then discuss the correction to the upper bound of the Lyapunov exponent caused 
by the non-uniformity of the PDF. In the following discussion, we assume that the intensity of external impulses is 
sufficiently small, and that the phase map F(8, e) is a strictly monotonically increasing function of 9, i.e. F'(6, e) > 0. 

A. Perturbative solution to the generalized Frobenius-Perron equation 

As a first step, we calculate the deviation of the stationary PDF from the uniform density perturbatively for small 
external impulses (up to the second order). Our starting point is the generalized Frobenius-Perron equation for the 
stationary PDF p(9), 

p{9)= f d(f>W{9 - 4>) f deQ{e) f di/>6(<f> - F(ip,e))p{tlj). (28) 
Jo J Jo 

We assume that the deviation of F(9,e) from the identity map F(9,e) = 9 is small, 

F(9,e) =9 + eG(9,e), (29) 

where we introduced a small parameter e in order to control the magnitude of the perturbation. By using a well-known 
formula for the 6- function, Ea. (|28fl can be rewritten as 

m - £ w - « / *g W 1F ^ e f e)l ■ (30) 

where ip*((f>, e) is a solution to F(ip, e) = <fi. We here used the fact that there exists only one solution, because F(if>, e) 
is a monotonically increasing function of ip (we do not consider the trivial case of F(tp, e) = ip, where p(9) always 
becomes uniform). 

We first calculate the solution ip*(4>,e) to F(ip,e) = 4> as a power series in e. We assume that the solution can be 
expanded in terms of e around the trivial solution tp* (</>, e) = <f> at e = as 

r (0, e) - + eVi (4>, e) + e V 2 * (<t>, e) + 0(e 3 ) . (31) 
By inserting this expression to F(tp*(<f>, e), e) = <j), we obtain 

<j> = 4> at 0(e°), 
V>i*(<?f>,e) + G(0,e) = at C^e 1 ), 
V> 2 *W>,e)+G'(0,eMW>,e) = at Q(e 2 ), 
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(32) 

Thus, to the second order in e, the solution ^}*(4>,e) is approximated by 

?/>*(0, e) = <f> - eG(<f>, e) + e 2 G'(0, e)G(</>, e) + 0(e 3 ). (33) 
Since F'(6, e) = 1 + eG'(9, e), it can be expanded as 

F'{ip*(cf>,e),e) = l + eG'(^e) - e 2 G"(0, e)G(0, e) + 0(e 3 ). (34) 
We also expand the stationary PDF p(9) in a power series of e as 

p(0) - 1 + e Pl {9) + e 2 P2 (6) + 0(e 3 ). (35) 

Since p(9) is normalized to 1, d9p\{9) = f d9p 2 {9) = should hold. Inserting the above expansions into Ea. (|30ll . 
we obtain 

l + ep 1 (9) + e 2 p 2 (9) + 0(e 3 ) 

= j d<t>W{e-4>) J deQ(e)[l + e{p 1 ( ( f > )-G'( ( f > ,e)} + e 2 {p 2 ( ( f,)~H^,e)} + 0(e 3 )] (36) 

where we utilized the fact that F'(ip*((f>, e), e) > 0, and defined 

H{4>, e) = P i(0)G(0, e) - G'(0, e)G(0, e). (37) 
Thus, the first order correction pi(9) to the uniform stationary PDF po(9) = 1 satisfies 

Pi(0) = ! d^W(d-cb) f deQ(e){p 1 (cf > )^G'(cl > ,e)} 



d<t>W(d-<P){ Pl (<P)-(G%(ct>)}, (38) 



and the second order correction p 2 (9) satisfies 



P2(0) 



f d<pW{9~<p) J deQ(e){p 2 (^)~H'(cj ) ,e)} 



d4W(e-<t>){p2(4>)-{H') e (<l>)}. (39) 



1 1 



Here and hereafter, for notational simplicity, we define an averaged function (f) e (0) of a function f(6,e) over Q(e) 
as 

(/> e (0) = / deQ(e)f(9,e), (40) 

such as (G') e (0) = / deQ(e)G' (4>, e) and (H') e ((f>) — J deQ(e)H'(<fi, e). We also define the Fourier transform between 
a function f{9) and its coefficient f(k) by 



oo „i 

f(9) = ^ ik6 f(k), f(k) = / , 

, JO 

K— — 00 



d9e- Z7!lkH f{9). (41) 

For example, the Fourier coefficients of pi(9) and G'(9, e) are denoted as pi(k) and G'(fc, e), respectively. The averaged 
Fourier coefficient (f) e (k) of f(k, e) over Q(e) is similarly defined as 

(/>e(fc)= / deQ(e)f(k,e), (42) 
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such as (G')e(m) = J deQ(e)G'(m, e) and (H') e {m) = f deQ(e)H'(m,e). 

Equations (|38|) and 139|) can be solved for px(9) and P2(#) through the Fourier transform, which yields 



px{m) = J deQ{e)W(m){p~i(m)-G'(m,e)} = W(m) {p~x(m) - (G') e (m)} , (43) 

and 

deQ(e)W(m) (p 2 {m) - H'(m,e)\ = W{m) |p 2 (m) - (ff') e (m)} . (44) 



p 2 {m) = 

It can easily be shown that the equations at m = give trivial relations, which should be neglected. We thus obtain 

W(m) ~ W(m) 
Pi m = \ j (GOe(m), p 2 (m) = \ J (ff')e(m) (45) 
Vy (m) — 1 Pv (m) — 1 

for m ^ 0, and the corrections pi(#) and p 2 {9) to the uniform density can be obtained as 

*W = E \ (gOe(m)e 2 ^, P2 (0) = 2 (g%(m)e 2 ^. (46) 

^( m ) - 1 ^( m ) - 1 

For the Poissonian impulses, the transition probability W{9) is given by Eq. Q21[l. and its Fourier coefficient IU(m) 
is given by 

r 1 a 

W(m) = / ddW(6)e- 27Time = (47) 

Jo A + 2irim, 

for integer m. Therefore, the coefficient in Ea. (|46fl is calculated as W(m)/[W(m) — 1] = — A/(2?rim). Using this, we 
can calculate the first order correction to the phase PDF as 

= -A[(G) e (tf)-G ], (48) 

where we defined G = (G) e (0) = /„ d9{G) e (9), and utilized the relation (G%(m) = (2 ■wim) (G) e (m) . Similarly, the 
second order correction to the PDF can be calculated as 

P^°) = E ^L^%(m) e 2 ™" 6 = -A£ We(m) e 2 " me 

= -A [(ff) e (0) - (H) e (0) , 

= -A [<pi(0)G(^», e)) e - <G'G) e (0)] + const. 

= A 2 [<G) e (#)] 2 - A 2 G • (G) e (0) + A(G'G) e (0) + const., (49) 

where we defined (G'G) e {9) = (G'(9, e)G(9, e)) e . The constant can be determined from the condition d9p 2 {9) = 0. 

Thus, in the Poissonian case, the averaged phase map (G) e (9) directly appears at the first order perturbation to 
the PDF, pi(9). Since A = 1/wt, the amplitude of the first order perturbation scales as e/ur, namely, the ratio of the 
impulse intensity to the non-dimensional time scale determined by the period of the limit cycle and the inter-impulse 
intervals. The second order perturbation p 2 (9) gives lowest-order nonlinear contributions from the phase map. 

B. Upper bound of the Lyapunov exponent 



The Lyapunov exponent A is bounded from above as 

A = ((logF'(^e))) e , e 
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d9p(9) J deQ(e)logF'(6,e) 
< J de P {0)lo g y deQ(e)F'(6,e) 



= J d6p{6) log |l + J deQ{e)eG'(9,e) 
< ddp(6) J deQ(e){eG'{9,e)} . 
Now, using Eq. (|35(l , correction to the upper bound of the Lyapunov exponent can be expanded as 
J d6p{8) J deQ(e){eG'{9,e)} 



(50) 



d9 {1 + e Pl (9) + ep 2 (8) + (9(e 3 )} / deQ(e) {eG'(8, e)} 



d9(G') e (9) + e 2 



d9p 2 (8)(G') e (9) + 0(e 4 ). 



(51) 



The first term corresponds to the (zero-th order) contribution from the uniform component of the phase PDF, which 
vanishes (similarly to the uniform case) irrespective of the functional form of the transition kernel W(9), 



d0(G') e (0) = (G%(1) - (G%(0) - 0. (52) 
For the Poissonian impulses, the first order correction to the upper bound of A can be calculated using Eq.ljJBJl as 

d9 Pl (6)(G') e (8) 
A f d0(G) e (G') e (6)+AG o I d6(G') e (6) 



-A 



{(G) e (0)} 2 



AG [{G) e (6)]l 



= 0. 



(53) 



Thus, the upper bound of the Lyapunov exponent A is still zero even if the first order correction to the uniform PDF 
is incorporated. The second order correction to the upper bound can similarly be calculated using Eg. 149(1 as 



d9p 2 (9){G') e (e) 



= A 2 



d9{(G) e (9)} 2 (G') e (9) - A 2 G / d9{(G) e (9)} (G%(6) 



+A d9(G'G) e (9)(G') e (9)+ const. d9{G') e {9) 



A 



2 r 



{(G%(9)Y 



-const. 



1 _ A 2 G r 
o 

i 



{{G')e(9)Y 



{(G') e (9)Y +A d9(G'G) e (9)(G') e (9) 



A / d9{G'G) e (6){G%(9). 



(54) 



Thus, only the term containing (G'G) e (9) in p 2 {9) gives non- vanishing contribution. Its sign cannot be determined 
at this point unless the explicit functional form of G(#, e) is given. This term could make the upper bound of the 
Lyapunov exponent slightly different from zero, but its effect is only of the order of e 3 . 
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Summarizing, the first order correction to the upper bound of the Lyapunov exponent A by the non-uniformity of 
the phase PDF is generally 0(e 2 ), but for the Poissonian impulses, it vanishes. Thus, the upper bound of A is still 
zero up to the first order approximation. The next order correction is only 0(e 3 ), which is quite small when e is small. 
Therefore, the impulse-induced phase synchronization will, in most cases, persist for weak Poissonian impulses even 
if the phase PDF becomes slightly non-uniform for small e. 



C. Examples of Lyapunov exponents 

1. Stuart- Landau oscillator 

Stationary phase PDFs p(8) of the Stuart-Landau oscillator driven by external impulses at a = 0.1 and r = 2 are 
shown in Figure El The curves are obtained by (i) a direct simulation of the original model Eq.Q without Gaussian- 
white noise, (ii) a direct simulation of the reduced phase model Eq. 1)1 5JI , (iii) a numerical calculation of the stationary 
solution of the corresponding Frobenius-Perron equation Ea.l|2ij|). and (iv) the perturbation theory up to the second 
order, respectively. All curves agree well, which confirms the validity of our approximation at least for small impulse 
intensity. In this case, the first order perturbation already gives a nice fit to the actual PDF, and the second order 
perturbation gives only a tiny correction. 

Figure plots the Lyapunov exponent A as a function of a at t — 2, which is obtained by (i) directly using 
the original model Eq.QJ without Gaussian- white noise, (ii) directly using the reduced phase model Ea. ((r3|l . (iii) a 
calculation using numerical phase maps and uniform phase PDF, and (iv) a calculation using numerical phase maps 
and the approximated phase PDF. Reflecting the symmetry of the limit cycle, the graph of A is also symmetrical 
with respect to a = 0. Since the phase map is always monotonically increasing in this range of <r, the Lyapunov 
exponent A calculated assuming uniform phase PDF (iii) is always non-positive and only becomes when a = 0. The 
Lyapunov exponent A calculated using approximate phase PDF (iv) is also always non-positive. Since the correction 
to the uniform PDF is small, the difference between (iii) and (iv) is also small. Both curves agree well with the actual 
Lyapunov exponent obtained by (i) and (ii). 

1.04 | . 1 1 1 1 1 . 1 1 1 



1.02 



1.00 



0.96 




• • Stuart-Landau 
x x Phase model 
— Frobenius-Perron 
-- Perturbation 



0.0 0.2 0.4 0.6 

e 



0.8 1.0 



FIG. 6: Stationary phase PDFs of the Stuart-Landau oscillator driven by external impulses obtained by (i) a direct simulation 
of the original model, (ii) a direct simulation of the corresponding phase model, (iii) numerically solving the corresponding 
Frobenius-Perron equation, and (iv) the perturbation theory. 



2. Hodgkin- Huxley model 

Similarly, stationary PDFs of the Hodgkin-Huxley neural oscillator driven by external impulses at a = 2 and 
t = 100 are shown in Figure|Sl As previous, the curves represent the results obtained by (i) a direct simulation of the 
original model, (ii) a direct simulation of the reduced phase model, (iii) a numerical solution of the Frobenius-Perron 
equation, and (iv) the perturbation theory. Of course, the results of (ii) and (iii) give a nice fit to the actual phase 
PDF obtained by (i). The result of perturbation theory (iv) also gives a reasonable fit to the actual phase PDF. As in 
the previous case, the first order perturbation gives a good fit to the actual PDF, and the second order perturbation 
gives only a tiny correction. 
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FIG. 7: Dependence of the Lyapunov exponent A on the impulse intensity a obtained by (i) directly using the original model 
Eq.@ without Gaussian- white noise, (ii) directly using the reduced phase model Eg. (1151 . (iii) a calculation using numerical 
phase maps assuming uniform phase PDF, and (iv) a calculation using numerical phase maps and the approximated phase 
PDF. 



Figure [5] plots the Lyapunov exponent A as a function of a at r = 100, which is obtained by (i) directly simulating 
the original model without external disturbances, (ii) directly simulating the reduced phase model, (iii) calculation 
using numerical phase maps assuming uniform phase PDF, and (iv) calculation using numerical phase maps and 
the phase PDF approximated up to the second order perturbation. Since the phase map is always monotonically 
increasing in this range of cr, the Lyapunov exponent A calculated assuming uniform phase PDF (iii) is always non- 
positive and only becomes when a = 0. The Lyapunov exponent A calculated using approximate phase PDF (iv) 
is also always non-positive. In this case, the correction to the uniform phase PDF is even smaller than the previous 
Stuart-Landau case, hence the Lyapunov exponents calculated by (iii) and by (iv) are almost indistinguishable. Of 
course, they coincide with the results obtained by direct simulation of the original model and the phase model. 

1.01 1 . 1 1 1 i 1 . 1 1 1 




FIG. 8: Stationary phase PDFs of the Hodgkin-Huxley neural oscillator driven by external impulses obtained by (i) direct 
simulation of the original model, (ii) direct simulation of the corresponding phase model, (iii) numerical solution of the cor- 
responding Frobenius-Perron equation, and (iv) perturbation theory using Fourier coefficients numerically obtained from the 
phase map. 



VI. SUMMARY 



We analyzed phase synchronization of general limit-cycle oscillators subject to external impulses by reducing the 
dynamics of the oscillator to a random phase map. We proved that when the phase maps are strictly monotonic and 
the mean inter-impulse interval of the input current is sufficiently large, the Lyapunov exponent of the system always 
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FIG. 9: Lyapunov exponent A obtained as a function of the impulse intensity a. by (i) directly using the original model, (ii) 
directly using the reduced phase model, (iii) a calculation using numerical phase maps assuming uniform phase PDF, and (iv) 
a calculation using numerical phase maps and the approximated phase PDF. 



becomes negative, leading to fluctuation-induced phase synchronization. We also treated the case where the inter- 
impulse interval is finite perturbatively for weak Poissonian impulses, and proved that the next order correction to 
the upper bound of the Lyapunov exponent is also zero, hence the fluctuation-induced phase synchronization persists 
even if the phase distribution becomes slightly non-uniform. 

Mathematically, the non-positivity of the Lyapunov exponent is a general result of the concavity of the log function 
and the monotonicity and periodicity of the phase map. Therefore, this result is not restricted to specific oscillators, 
but also holds generally for a wide variety of limit-cycle oscillators. Examining the significance of our results in 
practical problems would be an interesting topic. 

Though we did not derive in this paper, we can reduce the phase model driven by Poissonian impulses to an 
Ito-Langevin phase equation in the limit of weak and frequent impulses when the net drift induced by the external 
impulses vanishes. It yields 

6(t) = u} + Z(6)-ri(t), (55) 

where rj(t) is a n-dimensional Gaussian- white noise. However, it can be shown that the phase disturbance is neu- 
trally stable for this Ito-Langevin phase model, namely, the Lyapunov exponent is constantly zero |10|. as a direct 
consequence of the Ito stochastic integral j2^|. Therefore, if we take the Langevin limit within our phase model, 
fluctuation-induced phase synchronization does not occur. On the other hand, if the above Langevin phase equation 
is interpreted in the Stratonovich sense, which does not come out of the integration rule of the impulsive force we 
assumed in this paper, the phase synchronization occurs ^(j. Thus, slight difference in the treatment of the stochastic 
forcing leads to physically distinct results. Detailed discussions on this point, including the stochastic interpretation 
of impulsive forcing, will be reported in the future. 

Acknowledgments 

We thank D. Tanaka, J. Teramae, T. Aoyagi, and S. Nii for useful discussions. 



[1] Z. F. Mainen and T. J. Sejnowski, Science 268 1503 (1995). 

[2] R. R. de Ruyter van Steveninck, G. D. Lewen, S. P. Strong, R. Koberle, W. Bialek, Science, 275 1805 (1997). 

[3] Y. Tsubo, T. Kaneko. S. Shinomoto, Neural Networks 17 165 (2004). 

[4] R. V. Jensen, Phys. Rev. E. 58 R6907 (1998). 

[5] E. K. Kosmidis and K. Pakdaman, J. Comput. Neurosci. 14 5 (2003). 

[6] K. Pakdaman, Neural Comput. 14 781 (2002). 

[7] B. Gutkin, G. B. Ermentrout, and M. Rudolph, J. Comput. Neurosci. 15 91 (2003). 

[8] J. Ritt, Phys. Rev. E 68 041915 (2003). 

[9] J. M. Casado and J. P. Baltanas, Int. J. Bif. and Chaos 14 2061 (2004). 



18 



[10] J. N. Teramae and D. Tanaka, Phys. Rev. Lett. 93 204103 (2004). 
[11] K. Nagai, H. Nakao, and Y. Tsubo, Phys. Rev. E 71 036217 (2005). 

[12] D. S. Goldobin and A. Pikovsky, Phys. Rev. E 71 045201(R) (2005); Physica A 351 126 (2005). 

[13] A. T. Winfree, The Geometry of Biological Time (Springer- Verlag, New York, 2001, 1980). 

[14] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence (Springer- Verlag, Berlin, 1984). 

[15] A. Pikovsky, M. Rosenblum, and J. Kurths, Synchronization (Cambridge University Press, Cambridge, 2001). 

[16] D. L. Snyder, Random Point Processes (John Wiley & Sons, New York, 1975). 

[17] P. Hanggi, Z. Physik B 31, 407 (1978); Z. Physik B 36, 271 (1980). 

[18] A. Lasota and M. C. Mackey, Probabilistic properties of deterministic systems (Cambridge University Press, New York, 
1985). 

[19] L. M. Pecora and T. L. Carroll, Phys. Rev. A 44 2374 (1991). 
[20] A. S. Pikovsky, Phys. Lett. A 165 33 (1992). 

[21] P. Khoury, M. A. Lieberman, and A. J. Lichtenberg, Phys. Rev. E. 54 3377 (1996). 

[22] R. Toral, C.R. Mirasso, E. Hernandez-Garcia, O. Piro, in Unsolved Problems on Noise and Fluctuations, UPoN'99, D. 

Abbot and L. Kiss, eds. AIP Conference Proceedings 511, 255 (2000). 
[23] C. W. Gardiner, Handbook of Stochastic Methods (Springer, Berlin, 1997). 
[24] C. Koch, Biophysics of Computation (Oxford University Press, Oxford, 1999). 



